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I.  INTRODUCTION 

Nonlinear  optics  is  often  described  in  the  frequency  domain  so  that  dispersive  effects  can 
be  treated  in  a  simple  way.  The  convolution  integrals  associated  with  nonlinear  effects  are 
simplified  by  assuming  that  the  radiation  spectrum  takes  the  form  of  a  discrete  set  of  narrow 
peaks.  In  the  case  of  few  cycle  pulses,  the  convolution  integrals  cannot  be  reduced,  and  a  time 
domain  model  becomes  attractive.  The  finite-difference-time-domain  (FDTD)  technique  is 
a  well  established  method  for  solving  the  exact  Maxwell  equations  in  a  bounded  or  periodic 
domain  [1],  It  may  be  used  in  connection  with  the  particlc-in-ccll  (PIC)  technique  to  self- 
consistently  solve  for  the  motions  of  a  large  number  of  charged  particles  in  an  electromagnetic 
held  [2],  This  approach  is  often  used  to  model  fully  nonlinear  laser-plasma  interactions  and 
beam-plasma  interactions.  In  this  work,  it  is  extended  to  account  for  nonlinear  laser-crystal 
and  beam-crystal  interactions.  This  is  accomplished  by  incorporating  a  model  for  bound 
particles  into  the  PIC  code  turboWAVE  [3].  The  model  is  fully  parallelized,  and  runs  in  up 
to  three  dimensions. 

The  PIC  aspect  of  turboWAVE  has  much  in  common  with  a  number  of  other  codes 
designed  to  model  laser-plasma  interactions  [4-10].  The  nonlinear  optics  aspect  is  related 
to  a  number  of  other  codes  designed  to  model  ultrashort  pulse  propagation  in  a  nonlinear 
medium  [11-16].  Both  types  of  codes  can  be  divided  into  those  that  describe  the  radiation 
by  means  of  a  complex  envelope,  and  those  that  are  fully  explicit,  i.e.,  those  that  resolve 
the  optical  time  scale.  TurboWAVE  supports  both  models,  but  in  this  work  only  the  fully 
explicit  model  is  used.  As  a  result,  there  is  no  assumption  about  the  frequency  content  of  the 
radiation,  and  given  the  right  model  for  the  dielectric  response,  all  nonlinear  and  dispersive 
effects  are  accounted  for. 

Our  model  for  the  dielectric  response  generalizes  the  auxiliary  equation  technique  de¬ 
scribed  in  Refs.  [13,  15].  In  particular,  bound  charges  in  the  dielectric  are  represented  by 
effective  particles,  whose  contribution  to  the  four-current  is  computed  using  PIC  techniques. 
The  effective  particles  are  subjected  to  forces  arising  from  a  superposition  of  macroscopic  and 
microscopic  fields.  The  macroscopic  held  is  the  usual  electromagnetic  held  that  is  computed 
in  any  FDTD  code,  while  the  microscopic  held  is  a  three  dimensional  electrostatic  potential 
representing,  e.g.,  an  atomic  binding  potential.  The  resulting  equation  of  motion  is  expanded 
as  an  anharmonic  oscillator  equation.  By  using  multiple  species  of  effective  particles,  each 
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satisfying  a  different  oscillator  equation,  the  material  response  can  be  tailored  to  match  that 
of  real  materials  over  a  broad  range  of  frequencies.  Free  charges  are  self-consistently  incor¬ 
porated  into  the  model  by  linear  superposition  of  the  free  and  bound  sources  in  the  FDTD 
held  solver  [36].  The  resulting  model  is  capable  of  modeling  interactions  among  laser  pulses, 
particle  beams,  plasmas,  and  dielectric  crystals,  in  a  fully  dispersive,  nonlinear,  anisotropic, 
and  kinetic  way. 

In  Ref.  [17],  the  turboWAVE  extensions  described  in  this  report  were  introduced  for  the 
first  time.  In  Ref.  [18],  the  model  was  applied  to  a  novel  electro-optic  diagnostic  technique.  In 
the  present  report,  the  model  is  described  in  more  detail,  and  a  formal  analysis  of  numerical 
stability  is  given.  The  model  is  extended  to  include  third  order  nonlinearities,  and  the 
implementation  is  demonstrated  via  three  numerical  experiments. 


II.  DESCRIPTION  OF  THE  MODEL 


The  numerical  model  described  here  solves  the  exact  Maxwell  equations: 


r9D 

V  x  H  —  J  +  — — 

(la) 

<9B 

VXE  =  -^ 

(lb) 

VD  =  p 

(lc) 

V-B  =  0 

(Id) 

Here,  B  =  /U0H  +  M  and  D  =  eoE  +  P.  where  P  is  the  polarization  and  M  is  the  magneti¬ 
zation.  In  the  present  work,  we  take  M  =  0.  The  polarization  can  be  expressed  in  terms  of 
a  spatially  smoothed  four-current  due  to  bound  charges,  denoted  by  (( r /) ,  (j)).  This  results 
in  [19] 


1  (9E 

V  x  B  =  /i0  (J  +  (j))  + 

„  ^  <9B 

VxE=~m 

V  •  E  =  (p+  <»?))  /e0 


V  ■  B  =  0 


(2a) 

(2b) 

(2c) 

(2d) 
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With  this  formulation,  any  field  solver  designed  to  account  for  free  charges  can  be  easily 
adapted  to  account  for  bound  charges  by  making  the  substitution 


(p,  J)  -»■  (p+  (v) ,  J  +  (j» 


(3) 


Our  model  for  ((rj) ,  (j))  is  based  on  calculating  the  motions  of  effective  particles  that 
respond  to  the  superposition  of  a  microscopic  binding  potential  and  a  macroscopic  electric 
field  E.  Denoting  the  displacement  of  an  effective  particle  from  its  equilibrium  position  by 
r,  and  expanding  the  microscopic  potential  in  a  Taylor  series,  results  in  the  effective  particle 
equation  of  motion 
d2r 


dt 2 


jk  L 


dv ■ 

21 +  aijkrjrk  ~  br i'r yr j 


=  ~Ei 
m 


(4) 


where  the  subscripts  vary  over  Cartesian  coordinates,  and  q/m  is  the  charge  to  mass  ratio. 
The  anisotropy  of  the  medium  is  represented  by  the  tensors  T,  fi2,  and  a.  These  are, 
respectively,  the  damping  rate,  the  square  of  the  resonant  frequency,  and  the  first  anharmonic 
coefficient.  At  present,  the  second  anharmonic  coefficient,  b ,  is  assumed  scalar.  By  means 
of  a  coordinate  transformation,  Q2  can  always  be  made  diagonal.  The  basis  vectors  which 
induce  this  transformation  will  be  denoted  by  (eu,  e„,ew).  The  basis  vectors  used  for  the 
general  calculation  will  be  denoted  (ex,ey,ez).  The  crystallographic  basis  vectors,  which 
generally  coincide  with  (eu,  ev,  ew),  are  denoted  ((100) ,  (010) ,  (001)). 

In  general,  the  parameters  in  Eq.  (4)  depend  not  only  on  tensor  indices,  but  also  on  an 
index  p  that  identifies  a  particular  particle.  In  other  words,  each  particle  may  move  in  a 
unique  potential  well,  and  may  therefore  satisfy  a  unique  equation  of  motion.  The  system 
of  the  particle  together  with  its  equation  of  motion  is  called  an  oscillator.  In  practical 
implementations,  it  is  convenient  to  group  all  particles  that  satisfy  the  same  equation  of 
motion  into  an  oscillator  species  with  index  s.  Then  the  material  parameters,  q,  m,  Q2,  f~, 
a,  and  6,  depend  on  s,  while  only  r  is  explicitly  dependent  on  p.  As  discussed  below,  the 
macroscopic  sources  ((r/)  ,  (j))  can  be  derived  from  the  set  rp  using  a  variation  on  source 
deposition  techniques  developed  for  PIC  codes. 


III.  ANALYTICAL  DISPERSION  RELATION 

In  order  to  analyze  the  model  described  above,  consider  the  locally  averaged  displacement 
(rs)  associated  with  each  oscillator  species.  In  the  case  of  a  uniform  medium,  this  is  related 
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to  the  spatially  smoothed  four-current  by 


((h) ,  (j»  =  X]  V*N*  '  f  (r")  ’  (r«>) 


(5) 


where  Ns  is  the  density  of  oscillators,  and  the  tensor  f  is  an  oscillator  strength.  Assuming 
local  phase  coherence,  (rs)  satisfies  the  same  equation  as  rp,  which  in  the  (eu,ev,ew)  basis 

is 

w  +  2T4 + a‘‘-)  {r“]  =  ~^,{4  +  w<)'  (6) 

where  i  is  the  coordinate  index,  is  the  ith  coordinate,  A  is  the  vector  potential,  and  <f>  is 
the  scalar  potential.  Using  the  Coulomb  gauge,  the  radiation  is  described  by 


d 2 


qsNs  d  (ris)  d 24> 


+ 


e0  dt  did!;. 


(7) 


Consider  any  mode  with  a  wave-vector  that  is  collinear  with  one  of  the  principal  axes, 
(eu,e„,ew).  Then,  it  ce 
the  dispersion  relation 


(eu,e„,ew).  Then,  it  can  be  shown  that  (f>  —  0,  and  the  transverse  components  of  A  satisfy 


u 


p  2 

1  ,  Jiis^ps 


-  c2k 2  =  0 


(8) 


where  uj2  =  Nsqg/e0ms,  and 


T>ix(u)  =  Qjifl  —  2iuTns  —  ur 


(9) 


In  the  case  of  a  medium  comprised  of  a  single,  lossless,  isotropic  species  with  unit  oscillator 
strength,  the  dispersion  relation  is 

(O2  —  cu2)(c<;2  —  c2k2)  +  to2u2  =  0  (10) 

The  primary  feature  of  this  dispersion  relation  is  a  stop-band  between  the  resonant  frequency, 
O,  and  the  cutoff  frequency,  yTI2  +  u;2. 


IV.  OSCILLATOR  PARAMETERS 

Measured  dielectric  properties  are  usually  given  as  susceptibilities  in  the  frequency  do¬ 
main.  Hence,  in  order  to  apply  the  time  dependent  model  described  above  to  real  mate¬ 
rials,  it  is  necessary  to  connect  the  constants  fi2,  T,  f,  a,  and  b,  with  the  susceptibilities 
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Parameter 

Lattice  Oscillator  Electronic  Oscillator 

Oscillator  Strength,  / 

1 

1 

Resonance  Frequency,  R 

6.90  x  1013  rad/s 

6.38  x  1015  rad/s 

Damping  Frequency,  T 

6.25  x  1010  rad/s 

0 

Plasma  Frequency,  u jp 

9.27  x  1013  rad/s 

1.78  x  1016  rad/s 

Anharmonic  Coefficient,  0123 

-3.28  x  1037  m-1s-2 

4.1  x  1041  m”1s-2 

TABLE  I:  Oscillator  Parameters  for  Gallium  Phosphide 


u>2,  •••)•  In  the  (eu,  ev,  ew)  basis,  the  linear  susceptibility  has  only  the  diagonal  ele¬ 


ments 


p  2 

/  \  \  ^  J  iis^ps 

As 


Following  Ref.  [20],  the  first  anharmonic  coefficient  is  related  to  by 


1,U2,U3)  =  Y1 


ms  Vs  (oq )  Vs  (u2 )  Vs  (o>3 ) 


( aijk)s 


where  for  simplicity,  a  cubic  crystal  is  assumed.  The  nonlinear  susceptibility  is  related  to 
the  electro-optic  coefficient  Ty,*,  via 

+  Su,cj,  5uj)  =  -^  OzMOmfc(^,  5u)emj(u)  (13) 

lm 

Here,  the  factor  of  1/2  accounts  for  the  degeneracy  of  sum  and  difference  generation  when 
5u  — >  0,  and  e  is  the  relative  permittivity.  Finally,  the  nonlinear  refractive  index  is  given  by 


/  n  _  3  7/q  \  -  q2s _ UpS 

4  no(u>)2  m2s  U*(u>)Ds(u>)3 

where  no  is  the  linear  refractive  index  and  rj o  is  the  impedance  of  free  space  (see  appendix). 

Note  that  the  charge  to  mass  ratio  qs/ms  amounts  to  an  extraneous  free  parameter  that 

could  have  been  absorbed  into  qsNs ,  a.s,  and  bs.  In  this  work,  the  electronic  charge  to  mass 


ratio  is  always  assumed. 

As  an  example,  consider  the  cubic  crystal  gallium  phosphide  (GaP),  which  is  useful  for 
electro-optic  sensing  applications  [24-28].  Expressions  fitting  the  measured  linear  dispersion 
have  already  been  given  in  Ref.  [21].  There,  two  separate  expressions  were  used  for  the 
optical  and  THz  regimes.  A  two-oscillator  model  based  on  the  parameters  given  in  Table  I, 
agrees  with  both  expressions,  as  illustrated  in  Figs.  1(a)  and  (b). 


1*1231  (Pm/V) 
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FIG.  1:  Nonlinear  Lorentz  model  for  GaP.  Panel  (a)  displays  3i(e),  where  the  solid  curve  is  the 
two-oscillator  (20)  model  and  the  points  are  from  the  expressions  in  Ref.  [21].  Panel  (b)  is  similar 
to  (a)  except  9(e)  is  displayed.  Panel  (c)  displays  X123 (wl  ±  5ui,  u>l,  Su),  where  ujl  is  a  probe  laser 
frequency,  given  by  2i:c/ujl  =  0.6328  /xrn.  The  solid  curve  is  the  20  model,  and  the  dashed  curve 
is  the  expression  in  Ref.  [22] .  The  triangles  indicate  the  infra-red  frequencies  where  the  data  points 
in  [22]  are  located.  Panel  (d)  displays  r\23(uj,  5oo)  =  r4i(u> ,  5u) ,  where  5u>  — >  0.  The  curve  is  the 
20  model,  and  the  points  are  the  data  from  Ref.  [23].  The  20  model  cannot  be  made  consistent 
with  both  Refs.  [22]  and  [23]. 

Once  the  linear  dispersion  is  determined,  the  only  remaining  free  parameters  are  the 
anharmonic  coefficients  associated  with  each  oscillator.  In  order  to  choose  these  parameters, 
an  attempt  was  made  to  simultaneously  fit  data  from  Refs.  [22]  and  [23].  In  Ref.  [22], 
Xifk(uL  ±  Suj,  (jJl,  M  is  reported,  where  u>l  is  the  frequency  of  a  probing  helium-neon  laser, 
and  5u  is  one  of  several  infra-red  frequencies  produced  by  a  multi-line  H2+O2  laser.  In 
Ref.  [23],  rijk(co,6co)  is  measured,  where  5u  is  in  the  radio  frequency  (RF)  range,  and  u> 
takes  several  values  in  the  optical  range.  It  was  found  that  the  two-oscillator  model  cannot 
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be  made  to  agree  with  both  experiments.  The  parameters  of  Table  I  are  chosen  to  obtain 
agreement  with  the  measured  5u).  The  resulting  dispersion  is  displayed  in  Figs.  1(c) 

and  (d).  In  order  to  make  the  two-oscillator  model  agree  with  ±  8u,  8u),  a  much 

larger  absolute  value  of  the  lattice  anharmonic  coefficient  must  be  used.  This  leads,  roughly, 
to  the  oscillator  parameters  given  in  Ref.  [17],  although  for  the  simulations  presented  there 
and  in  Ref.  [18],  the  lattice  oscillator  was  taken  to  be  linear. 


V.  NUMERICAL  TECHNIQUE 


A.  Oscillator  Equation  and  Source  Deposition 


Most  electromagnetic  PIC  codes  use  the  leap-frog  technique.  More  specifically,  E'"+1  is 
computed  using  En,  Bn+2  and  where  n  is  the  time  level.  Then,  Bn+3  is  computed 

using  Bn+  2  and  En+1.  In  some  cases  a  Poisson  solver  is  used  to  refine  En+1  using  pn+l .  As 
discussed  above,  any  such  held  solver  can  be  used  in  the  presence  of  an  oscillator  population 
by  making  the  substitution  (p,  J)  — > ►  (p  +  (p)  ,  J  +  (j)).  This  implies  that  (p)  is  known  at 
integral  time  levels,  and  (j)  is  known  at  half-integral  time  levels. 

To  compute  (p)n+1  and  (j ) n+  2 ,  the  displacements  r”  and  r™+1  are  needed.  Let  T  repre¬ 
sent  the  matrix  of  transformation  of  a  vector  from  the  (eu,ev,ew)  basis  to  the  (ex,ey,ez) 
basis.  If  all  quantities  characterizing  the  oscillators  are  stored  in  the  (eu,  e,, ,  ew)  basis,  the 
displacements  can  be  updated  using 


^n+1 


2  rf  +  r. 


n—  1 


At2 


+  2  Tj 


rn+l 


Kn 


At 


+  (fi  V"  =  ^ 


m 


(15) 


where  At  is  the  time  step  and 


F"  =  qT^E"  -  marnrn  +  mb{ rn  •  rn)rn 


(16) 


The  macroscopic  four-current  ((p) ,  (j))  associated  with  the  microscopic  orbits  r”  can  be 
computed  using  PIC  techniques.  These  techniques  regard  each  particle  as  a  “cloud”  with 
dimensions  on  the  order  of  a  cell  size.  Knowing  the  orbit  of  each  cloud’s  centroid  allows 
one  to  compute  the  amount  of  charge  in  each  cell  as  a  function  of  time,  as  well  as  the 
current  passing  through  each  cell  wall.  Various  schemes  have  been  developed  to  facilitate 
this  calculation  [2,  29,  30].  In  the  case  of  bound  charges,  these  schemes  have  to  be  applied 
with  care  in  order  to  minimize  round-off  errors  associated  with  very  small  displacements. 


FIG.  2:  Charge  Deposition,  (a)  Geometry  of  a  grid  cell  in  relation  to  the  points  on  the  Yee 
mesh.  The  three  dimensional  arrows  point  to  the  mesh  points  where  the  corresponding  electric 
field  components  are  known,  (b)  Quadratic  Weighting.  The  bound  charges  are  represented  by 
charge  clouds  (solid  curve)  whose  equilibrium  position  (dashed  curve)  is  the  centroid  of  the  cell 
volume.  For  an  electron,  the  charge  in  a  cell  is  increased  (decreased)  by  the  area  of  the  intersection 
of  the  blue  (red)  area  with  the  cell. 


To  deposit  the  sources  due  to  bound  charges,  place  one  oscillator,  per  species,  in  each 
grid  cell  [37].  The  cell  geometry  is  shown  in  Fig.  2(a).  The  deposition  of  charge,  in  one 
dimension,  is  illustrated  in  Fig.  2(b).  The  oscillator  is  represented  by  two  oppositely  charged 
density  distributions.  The  equilibrium  distribution,  shown  as  a  dashed  line,  is  immobile. 
The  displaced  distribution,  shown  as  a  solid  line,  executes  the  orbit  described  by  Eq.  (15). 
The  net  charge  in  a  cell  is  the  integral  over  the  cell  volume  of  the  difference  between  the 
distributions.  To  minimize  round-off  error,  this  integral  should  not  be  carried  out  by  re¬ 
using  a  PIC  routine  that  deposits  each  distribution  separately.  Instead,  expressions  should 
be  written  for  the  net  charge  in  a  cell,  and  expanded  to  lowest  order  in  the  displacement. 
Then,  for  the  oscillator  in  cell  ( i,j ,  k),  the  net  charge  density  is  distributed  in  the  surrounding 
cells  as  follows: 


= 

(v)lj±1,k  = 
fa>u,*±l  = 


I  -rr  n-^sQs 

±er  ■  Tfr"— — 

x  p2Ax 

±e  •  Tfrn— ^ 

y  p2Ay 


±e,  •  Tfr1 


t NsQs 

'  2Az 


(17a) 

(17b) 

(17c) 


Here,  the  cell  dimension  is  Ax  x  Ay  x  Az.  The  current  density  through  the  surrounding 
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cell  walls  follows  directly  from  charge  conservation: 


“b  ,j,k 


(j) 


,,TJ 

i,j,k±  5 


±  (<C±+A  -  »W)  ^ 

±  -  «*«)  If 

±  (W$±i  -  Whfcti)  |f 


(18a) 

(18b) 

(18c) 


Note  that  the  above  expressions  are  particular  to  the  choice  of  quadratic  weighting.  This 
choice  provides  adequate  smoothing,  while  also  allowing  for  abrupt  vacuum- dielectric  tran¬ 
sitions.  A  useful  fact  is  that  the  current  density  due  to  an  effective  particle,  through  any 
cell  wall  enclosing  that  particle,  is  the  normal  component  of  | Nsqsvp ,  where 


vp  =  Tf 


„n+l  _  _n 

_P _ _P_ 

At 


(19) 


is  the  velocity  of  the  effective  particle  in  the  (ex,ey,ez)  basis. 


B.  Linear  Stability  and  Accuracy 

In  an  ordinary  PIC  code,  the  primary  stability  criterion  is  the  Courant  condition,  which 
in  one  dimension  reads  cAt  <  A z,  where  At  is  the  time  step  and  Az  is  the  grid  spacing. 
It  is  clear  that  in  the  presence  of  a  dispersionless  dielectric,  this  would  be  modified  to 
read  cAt/n  <  A z,  where  n  is  the  refractive  index.  In  the  case  of  a  fully  dispersive  model 
such  as  the  one  described  above,  the  modification  of  the  Courant  condition  is  more  subtle. 
Expressing  Eqs.  (6)  and  (7)  in  finite  difference  form,  assuming  propagation  along  a  principal 
axis,  and  inserting  exponential  forms  for  Ai  and  (r*s),  gives  the  numerical  dispersion  relation 


where 


and 


cos  u  At  —  V(k)  +  Vi(  oS) 


V(k)  =  1  -  2c2 At 


n») = E  t, 


2a^2  i  sin2  7,kx Ax  sin2  \kyAy  sin2 1 kzAz 


+ 


+ 


Ax2  Ay 2 

sin2  t 


A  z* 


(20) 


(21) 


(22) 


A/-  -  (1  +  riisAt)  sin2  l^A t  -  (,  I  ...  A/  -m 
Note  that  when  the  density  of  oscillators  vanishes,  Vi  — >  0  and  the  well  known  vacuum 
equation  is  recovered.  The  effect  of  Vi  can  be  estimated  by  assuming  Vt2At2  -C  1  and 
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kx  Ax  rs_/  kyAy  rs_/  k%Az  ^  c o  At  ^  7 r.  These  assumptions  are  1 1 lotiva ted  by  tliG  fa c t  I  lull 
the  resonance  frequency  should  be  well  resolved,  and  by  the  expectation  that  numerical 
instabilities  are  most  severe  at  the  Nyquist  frequency.  Under  these  assumptions,  the  stability 
criterion  becomes 


At  < 


c2  c2  c2  y1/2 
Arc2  +  Ay2  +  A22 ) 


(1  _  1  fas^psAt2  \  1/2 

y  4  1  +  viisAt  J 


(23) 


This  is  more  stringent  than  the  usual  Courant  condition.  Numerical  experiments  show  that 
Eq.  (23)  is  accurate  in  practice. 

In  order  to  evaluate  the  accuracy  and  stability  of  the  differencing  scheme  more  precisely, 
the  numerical  dispersion  relation  is  evaluated  numerically  for  the  case  of  a  single,  lossless 
oscillator  species.  Since  the  dispersion  relation  scales,  we  normalize  all  frequencies  to  the 
resonant  frequency,  12,  taking  up  =  212,  Az  =  0.5c/12,  and  Arc  =  Ay  — >  00.  The  results 
are  shown  in  Fig.  3  for  three  choices  of  the  time  step.  The  case  At  =  0.2/12  is  shown  in 
panel  (a),  the  case  At  =  0.4/12  is  shown  in  panel  (b),  and  the  case  At  =  0.48/12  is  shown 
in  panel  (c).  The  analytical  dispersion  relation  is  overlayed  for  comparison.  As  usual,  the 
best  accuracy  is  obtained  for  small  values  of  uj  and  k.  As  cAt/Az  increases,  the  accuracy 
improves,  until  numerical  instability  sets  in.  This  behavior  is  similar  to  the  vacuum  case. 
However,  unlike  the  vacuum  case,  instability  occurs  for  cAt/ Az  <  1,  as  shown  in  panel  (c). 
Better  accuracy  may  be  obtained  by  decreasing  both  Ax:  and  At  proportionately. 


C.  Nonlinear  Stability 


In  the  case  of  the  fully  nonlinear  system,  a  rigorous  stability  analysis  is  difficult.  However, 
if  the  oscillator  equation  is  viewed  as  describing  an  effective  particle  in  a  potential  well,  the 
stability  criterion  may  be  viewed  as  the  requirement  that  the  particle  never  enter  any  region 
of  runaway  acceleration.  Such  a  region  exists  if  b  >  0,  for  in  this  case  there  is  a  divergent 
repulsive  force  as  r  — »  ±00.  Instability  is  also  possible  if  b  =  0  and  a  /  0.  Assuming  there 
are  no  metastable  regions,  the  stability  criterion  can  be  written 

Fnl  •  r  <  mQ.2  r  •  r  (24) 

where 

Fnl  =  mb(r  ■  r)r  —  marr  (25) 
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FIG.  3:  Numerical  dispersion  for  (a)  cAt/Az  =  0.4,  (b)  cAt/Az  =  0.8,  and  (c)  cAt/Az  = 
0.96.  Solid  curves  are  the  analytical  dispersion  relation,  solid  circles  are  3?(u;/fi),  and  open  circles 
connected  by  a  thin  line  are  3f(u;/fi).  In  all  cases,  ujp/VI  =  2,  and  A z  =  0.5 c/Q.  Panel  (c)  shows 
that  numerical  instability  may  occur  for  cAt/Az  <  1. 


If  one  is  interested  primarily  in  second  order  effects,  the  system  can  be  stabilized  by  taking 
b  <  0.  The  magnitude  of  b  is  chosen  such  that  the  order  of  magnitude  of  the  first,  second, 
and  third  order  forces  are  the  same  at  some  critical  point.  This  leads  to 

,2 


,  a~ 
b  ss - 

d2 


(26) 


where  a  and  hi  are  typical  values  of  the  elements  of  a  and  fi.  If  third  order  effects  are 
important  physically,  a  stabilizing  fifth  order  nonlinearity  can  be  easily  introduced.  In  this 
case, 

Fnl  =  md( r  •  r)2r  +  mb( r  •  r)r  —  marr  (27) 


where  d  is  the  fifth  order  anharmonic  coefficient.  Demanding  that  the  first,  third,  and  fifth 
order  forces  be  equal  at  some  critical  point  gives 


d 


D2 


(28) 


The  ultimate  nonlinear  stabilization  technique  would  utilize  physical  effects,  such  as  ion¬ 
ization  in  the  case  of  the  electronic  oscillator,  or  material  damage  in  the  case  of  the  lattice 
oscillator.  However,  accounting  for  these  processes  is  beyond  the  scope  of  the  present  work. 


Re(€) 
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Parameter 

Lattice  Oscillator 

False  Resonance 

Oscillator  Strength,  / 

1 

1 

Resonance  Frequency,  If 

6.90  x  1013  rad/s 

8.91  x  1014  rad/s 

Damping  Frequency,  T 

6.25  x  1010  rad/s 

0 

Plasma  Frequency,  up 

9.27  x  1013  rad/s 

2.46  x  1015  rad/s 

Anharmonic  Coefficient,  0123 

-3.28  x  1037  m-1s-2 

1.23  x  1040  m-1s-2 

TABLE  II:  False  Resonance  Model  for  Gallium  Phosphide 


6w!2n  (THz) 


FIG.  4:  False  resonance  model  for  GaP.  Panel  (a)  displays  3ft(e),  where  the  solid  curve  is  the 
false  resonance  model  and  the  points  are  from  the  expressions  in  Ref.  [21].  Panel  (b)  displays 
<W,  Slu),  where  the  curve  is  the  false  resonance  model  and  the  points  are  the  true  resonance 
model  (Table  I). 


D.  False  Resonance  Technique  for  Accelerating  Computation 

It  sometimes  happens  that  the  accuracy  and  stability  criteria,  Q2At2  <C  1  and  uJ2At2  -C  1, 
force  one  to  use  a  very  large  number  of  time  steps  to  simulate  a  given  process.  This  limitation 
can  be  mitigated  if  the  frequency  range  of  interest  is  well  below  O.  In  such  cases,  it  is 
possible  to  artificially  shift  0  to  a  lower  frequency  while  preserving  the  dispersion  relation 
in  the  region  of  interest.  In  this  way,  the  “true  resonance”  is  replaced  by  a  “false  resonance” 
which  is  easier  to  resolve,  but  leads  to  the  same  physics. 

As  an  example  of  a  false  resonance  model,  consider  again  GaP,  which  has  the  true  reso¬ 
nance  parameters  displayed  in  Table  I.  When  modeling  an  electron  bunch  with  a  characteris¬ 
tic  time  of,  say,  100  femtoseconds,  the  frequency  range  of  interest  is  well  below  the  electronic 
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resonance,  which  has  2ti /Q,  ~  1  femtosecond.  One  may  then  use  a  false  resonance  model 
such  as  the  one  described  by  Table  II.  The  linear  and  nonlinear  dispersion  curves  associated 
with  this  model  are  shown  in  Fig.  4.  In  the  THz  range,  the  dispersion  relation  remains 
accurate.  Note  that  the  dispersion  in  this  range  could  be  fine  tuned  by  adding  additional 
false  resonances. 


E.  Boundary  Conditions  and  Moving  Window 

In  simulations  of  laser-plasma  interactions,  the  moving  window  technique  [31]  is  often 
used.  In  this  technique,  the  mesh  points  are  shifted  by  one  cell  in  between  time  levels 
such  that  the  computational  region,  or  “window,”  moves  at  the  speed  of  light.  Boundary 
conditions  are  greatly  simplified  by  the  requirements  of  causality.  When  the  plasma  density 
is  low,  the  laser  pulse  group  velocity  is  nearly  the  speed  of  light,  and  the  pulse  stays  in  the 
window  for  a  long  time.  Hence,  the  window  only  has  to  be  as  long  as  the  laser  pulse,  even 
if  the  interaction  region  is  much  longer. 

When  the  group  velocity  is  significantly  smaller  than  the  speed  of  light,  the  pulse  quickly 
falls  behind  a  light  frame  window,  and  one  might  as  well  work  in  the  lab  frame.  If  the 
window  speed  is  the  group  velocity,  a  fully  dispersive  model  might  allow  precursor  signals 
to  reach  the  front  of  the  window.  Moreover,  signals  reflected  from  the  back  of  the  window 
might  work  their  way  toward  the  region  of  interest.  One  solution  is  to  use  perfectly  matched 
layers  (PML)  [32]  to  absorb  waves  near  the  boundaries  of  the  window.  In  order  for  this  to 
be  effective,  it  is  necessary  to  use  a  large  number  of  layers  due  to  the  fact  that  the  moving 
window  involves  shifting  the  field  data  by  one  cell  between  time  levels.  In  order  that  this 
shift  should  weakly  perturb  the  solution,  the  PML  conductivity  has  to  change  very  gradually 
from  one  cell  to  the  next. 

In  the  case  of  lab  frame  simulations,  another  option  is  to  use  Lindman’s  absorbing  bound¬ 
ary  condition  [33]  in  the  longitudinal  direction,  and  periodic  boundary  conditions  trans¬ 
versely.  If  a  Poisson  equation  has  to  be  solved,  the  solution  in  the  box  can  be  matched  to  a 
decaying  solution  outside  the  box  [30]. 
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VI.  NUMERICAL  EXPERIMENTS 

The  model  described  above  has  been  implemented  as  a  module  within  the  turboWAVE 
[3]  framework.  In  this  section,  three  numerical  experiments  are  presented.  First,  the  electro- 
optic  effect  in  GaP  is  demonstrated  by  comparing  the  numerical  change  in  polarization  with 
that  predicted  by  the  low  frequency  theory.  Second,  soliton  propagation  in  fused  silica  is 
demonstrated.  Finally,  to  take  full  advantage  of  the  PIC  framework,  a  relativistic  electron 
bunch  is  passed  near  a  GaP  crystal,  and  the  induced  fields  are  examined. 

A.  Electro-Optic  Effect 

In  this  numerical  experiment,  an  x-polarized  laser  pulse  is  copropagated  with  a  ^/-polarized 
THz  half-wave  in  a  GaP  crystal,  and  the  change  in  the  polarization  state  of  the  laser  pulse  is 
measured.  The  length  of  the  THz  half-wave  is  chosen  to  be  long  enough  so  that  the  standard 
picture  of  the  electro-optic  effect  applies.  In  this  picture,  the  change  in  the  impermeability 
[38]  tensor  due  to  a  DC  electric  field  is 

A Tjij  r  ijkEk  (29) 

For  any  Zinc-Blende  type  crystal,  rVJk  =  r 123  =  r 41  if  all  three  indices  are  distinct,  and 
Tijk  =  0  otherwise.  The  crystal  orientation  is  represented  by  the  operator  T,  which  can  be 
viewed  as  a  sequence  of  rotations  that  re-orients  the  crystal  starting  with  the  (eu,  e„,  e<y) 
and  (ex,  ey,  ez)  bases  aligned.  It  turns  out  that  the  electro-optic  effect  is  maximized  if  [34] 

T  =  Rz(ir/2)Rx(-n/2)Rz(ir/4)  (30) 

with  the  laser  field  polarized  in  the  x-direction  and  the  DC  field  polarized  in  the  ^/-direction. 
Here,  Ri  is  an  active  right-handed  rotation  about  the  ith  coordinate  axis.  As  an  example, 
if  E  and  Ary  are  given  in  the  (ex,  ey,ez)  basis,  and  r  is  given  in  the  (eu,ev,ew)  basis,  then 
Ary  =  TrT_1ET_1.  To  calculate  the  change  in  polarization  due  to  the  electro-optic  effect,  a 
principal  axis  transformation  is  usually  made  [34],  Alternatively,  one  may  make  direct  use 
of  the  solution  of  the  slowly  varying  envelope  equation, 

£(z)  =  exp  £(z  =  0)  (31) 


15 


where  £  is  the  complex  electric  held  envelope,  Ao  is  the  free  space  laser  wavelength,  n  is  the 
refractive  index  of  the  unperturbed  crystal,  and  Ae  is  the  change  in  permittivity  induced  by 
the  DC  held.  This  equation  remains  valid  even  when  Ae  is  a  tensor  and  £  is  a  vector,  and 
is  easily  evaluated  in  a  software  environment  that  supports  the  exponential  of  a  matrix. 

We  now  compare  the  prediction  of  Eq.  (31)  with  the  results  of  a  turboWAVE  simulation. 
In  the  simulation,  an  x-polarized  laser  pulse  and  a  //-polarized  THz  half-wave  are  incident  on 
a  GaP  crystal  oriented  as  in  Eq.  (30).  The  simulation  parameters  are  given  in  Tables  I  and 
III.  The  lattice  nonlinearity  was  stabilized  with  b  =  — a2/£l 2 .  The  results  from  the  simulation 
are  shown  in  Fig.  5.  The  data  are  evaluated  60  /mi  into  the  crystal  to  avoid  interference  from 
reflections  generated  at  either  vacuum- crystal  interface.  Since  the  turboWAVE  model  uses 
a  universal  electric  held,  the  optical  and  THz  frequency  components  have  to  be  extracted 
by  means  of  a  Fourier  hlter.  The  envelope  of  the  optical  wave  is  constructed  by  applying  a 
hard-edged  band-pass  hlter  centered  at  the  optical  frequency,  +u;0.  The  result  is  then  shifted 
in  frequency  by  —  cn0,  multiplied  by  2  (to  account  for  the  energy  in  the  negative  frequencies) 
and  the  inverse  Fourier  transform  is  applied.  Based  on  Fig.  5,  the  polarization  ratio  at 
z  —  60  /an,  and  at  times  prior  to  the  generation  of  internal  rehections,  is  \£y/£x\2  ~  0.0018. 
Applying  the  analytical  formula  with  the  inputs  Ey  =  68  kV/cm,  n  =  3.16,  and  r'123  =  0.89 
pm/V  gives  \£y/£x\2  ~  0.0019.  Thus,  the  simulation  model  closely  agrees  with  the  theoretical 
prediction. 

An  interesting  feature  of  the  simulation  data  is  that  the  second  //-polarized  optical  pulse  is 
larger  than  the  hrst.  This  appears  to  happen  because  the  rehected  x-polarized  optical  pulse 
continues  to  be  rotated  by  the  reflected  //-polarized  THz  pulse.  This  leads  to  a  coherent 
reinforcement  of  the  reflected  //-polarized  optical  pulse. 

B.  Soliton  Propagation 

In  this  numerical  experiment,  a  20  femtosecond,  2.1  /mi  pulse  is  propagated  in  fused  silica 
in  both  the  linear  and  soliton  regimes.  The  Sellmcier  formula  for  fused  silica  can  be  matched 
exactly  to  a  three  oscillator  model,  the  parameters  of  which  are  displayed  in  Table  IV.  The 
second  anharmonic  coefficient  is  chosen  to  give  a  nonlinear  refractive  index  of  n2  =  3  x  10“ 16 
cm2/W.  The  choice  to  take  b  ^  0  only  in  the  third  oscillator  has  no  effect  other  than  on  the 
frequency  dependence  of  n2. 
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Parameter 

Symbol 

Value 

Time  Step 

cAt 

0.012  /xm 

Space  Step 

Az 

0.017  /xm 

Cells 

Nz 

213 

Steps 

Nt 

217 

THz  Field 

Ey 

68  kV / cm 

THz  Half  Wave 

th 

2.2  ps 

Laser  Field 

I  f  I 

18  kV/crn 

Laser  FWHM 

tl 

220  fs 

Laser  wavelength 

Ao 

0.81  /xm 

Crystal  Length 

L 

120  /xm 

Crystal  Orientation 

T 

Eq.  30 

TABLE  III:  Parameters  for  simulation  of  electro-optic  effect.  The  fields  are  measured  inside  the 
dielectric.  The  THz  field  is  a  half-cycle  pulse  with  base-to-base  width  th- 


Time  (ps) 


FIG.  5:  Simulated  electric  field  vs.  time  at  2  =  60  /xm  from  the  entrance  plane  of  the  crystal, 
(a)  //-polarized  THz  field  (low-pass  Fourier  filter)  (b)  x-polarized  optical  field  envelope  (band-pass 
Fourier  filter  with  shift)  (c)  //-polarized  optical  field  envelope  (band-pass  Fourier  filter  with  shift). 
The  incident  optical  field  is  a  single,  purely  x-polarized  pulse,  and  the  incident  THz  field  is  a  single, 
purely  //-polarized  half-wave.  The  additional  optical  pulses,  and  the  late-time  distortion  of  the  THz 
field,  are  due  to  reflections. 


17 


Parameter 

Osc. 

1 

Osc. 

2 

Osc.  3 

Unit 

Oscillator  Strength,  / 

1 

1 

1 

Resonance  Frequency,  f! 

2.35  x  10 

14 

1.63  x  10 

16 

2.78 

X 

1016 

rad/s 

Damping  Frequency,  T 

0 

0 

0 

rad/s 

Plasma  Frequency,  uop 

1.79  x  10 

14 

1.06  x  10 

16 

2.30 

X 

1016 

rad/s 

Anharmonic  Coefficient,  l 

) 

0 

0 

4.72 

X 

1054 

m~2s-2 

TABLE  IV:  Oscillator  Parameters  for  Fused  Silica 


Parameter 

Symbol 

Value 

Time  Step 

cAt 

0.014  pm 

Space  Step 

Az 

0.034  pm 

Cells 

Nz 

212 

Steps 

Nt 

4  x  105 

Window  Speed 

V9 

0.68c 

PML  thickness 

- 

8.7  pm 

Intensity 

To 

0.4  TW/crn2 

Pulse  Width 

To 

20  fs 

Free  Space  Wavelength 

Ao 

2.1  pm 

Dispersion  Length 

Ld 

2.8  mm 

Nonlinear  Length 

Lnl 

2.8  mm 

GVD  Coefficient 

02 

—  1.42  x  10~25  s2/m 

TABLE  V:  Parameters  for  simulation  of  soliton  propagation. 

The  amplitude  envelope  of  the  fundamental  soliton  has  the  form  sech(r),  where  r  = 
(t  —  z/vg)/T0,  vg  is  the  group  velocity,  and  T0  is  the  pulse  width  parameter.  In  simulations, 
it  is  convenient  to  have  a  function  whose  support  is  strictly  in  a  finite  region.  Furthermore, 
in  a  model  of  the  type  used  here,  the  initial  conditions  depend  on  z,  and  the  pulse  starts  in 
vacuum.  Hence,  the  appropriate  initial  condition  is  [39] 

Ex{t  =  0,  z)  =  e1/1|£’o|sech(C/cT0)  sm(u0(/c)C((/cT0)  (32) 


where  (  =  z  —  zq,  Zq  fixes  the  spatial  origin,  and  C(x)  is  a  suitable  clipping  function.  The 
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FIG.  6:  Simulation  of  ID  soliton  propagation  in  fused  silica.  Plots  show  e1^2E2/r]Q  vs.  vgt  —  z , 
where  vgt  is  held  fixed  in  each  plot,  and  vg  is  the  velocity  of  the  moving  window.  Panel  (a)  shows 
the  initial  pulse,  (b)  shows  the  pulse  after  8  mm  propagation,  and  (c)  shows  the  result  from  (b)  in 
the  case  where  the  initial  intensity  is  Jo/ 100. 


clipping  function  used  here  is 


C{x) 


\x\  <  4 
<  4  <  |x|  <  6 
Id  >  6 


1 

cos2[(|a;|  —  4)7t/4] 

0 


(33) 


The  parameters  of  the  pulse,  and  the  numerical  parameters  associated  with  the  simulation, 
are  displayed  in  Table  V.  The  intensity  and  electric  field  are  related  by  I0  =  e1/2|£0|2/2/7o- 
The  free  space  wavelength  is  A0  =  2ttc/uj0.  The  group  velocity  dispersion  (GVD)  coefficient 
is  defined  by  c/^2  =  d2(e1^2uj) /dt2]^^.  The  dispersion  length,  Ld  =  T02/ 1/32|,  measures  the 
propagation  length  needed  to  observe  group  velocity  dispersion  (GVD)  effects.  The  nonlinear 
length,  Lml  =  c/ujo^Io,  measures  the  propagation  length  needed  to  observe  nonlinear  pulse 
distortions.  The  conditions  j32  <  0  and  LD  =  Lnl  result  in  stable  propagation  of  the  pulse 
without  distortion  [20,  35].  This  is  illustrated  in  Fig.  6.  Comparison  of  panels  (a)  and  (b) 
shows  that  when  /0  =  0.4  TW/cm2,  the  pulse  maintains  its  original  form  after  propagating 
nearly  3 Lp.  Panel  (c)  shows  that  when  I0  is  reduced  by  a  factor  of  100,  the  pulse  broadens 
significantly  over  the  same  propagation  distance. 


19 


Parameter 

Symbol 

Value 

Time  Step 

cAt 

0.168  y m 

Space  Step 

< 

II 

a 

<1 

II 

<1 

0.672  ym 

Cells 

Nxx  Nyx  Nz 

to 

0 

X 

to 

0 

X 

to 

0 

Steps 

Nt 

3  x  104 

Window  Speed 

V 

c 

Bunch  Diameter 

n 

8  ym 

Bunch  Length 

n 

80  fs 

Electron  Energy 

(7  —  l)mc2 

250  MeV 

Crystal  Size 

Lx  x  Ly  x  Lz 

605  x  360  x  00  ym 3 

Crystal  Orientation 

T 

Eq.  30 

TABLE  VI:  Parameters  for  simulation  of  fields  induced  by  an  electron  bunch. 

C.  Electron  Bunch  Passing  Near  a  Crystal 

A  three  dimensional  simulation  of  the  fields  generated  in  an  electro-optic  crystal  by  a 
passing  relativistic  electron  bunch  was  described  in  Ref.  [17].  Here,  a  similar  simulation  is 
carried  out,  with  the  following  distinctions.  First,  the  bunch  charge  is  varied  until  significant 
nonlinear  distortions  are  observed  in  the  bunch  fields.  Second,  the  false  resonance  technique 
is  used  to  accelerate  the  calculation,  rather  than  the  sub-cycling  technique  that  was  used 
in  Ref.  [17].  The  former  technique  is  found  to  have  superior  stability  and  conservation 
properties.  The  false  resonance  oscillator  parameters  are  the  same  as  those  given  in  Table  II, 
except  that  a  stabilizing  third  order  nonlinearity  is  added  to  each  oscillator.  The  other 
parameters  of  the  simulation  are  given  in  Table  VI. 

Figs.  7(a)  and  (b)  display  the  electric  held  component  Ey(x  =  0 ,y,z),  for  the  cases 
Q  =  0.28  pC  and  Q  =  2.8  pC,  where  Q  is  the  bunch  charge.  Figs.  7(c)  and  (d)  show  the 
corresponding  line-outs  at  y  =  0.  The  bunch  charge  affects  the  degree  of  nonlinearity  in 
the  interaction.  The  primary  features  are  the  vertical  phase  fronts  located  in  the  range 
—400  <  z  —  ct  <  —300  /mi,  and  the  diagonal  phase  front  that  appears  at  all  z  —  ct.  As 
discussed  previously  [17],  these  can  be  roughly  associated  with  coherent  transition  radiation 
(CTR)  and  Cherenkov  radiation,  respectively.  In  electro-optic  sensing  applications,  the 
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FIG.  7:  Simulation  of  fields  induced  by  an  electron  bunch  with  (a,c)  Q  =  0.28  pC  and  (b,d)  Q  =  2.8 
pC.  The  three  dimensional  data  is  evaluated  at  x  =  0  in  order  to  present  a  two  dimensional  plot. 
The  green  oval  represents  the  location  of  the  electron  bunch,  which  propagates  from  left  to  right. 
The  boundary  between  the  crystal  and  vacuum  is  located  at  y  =  200  /im.  The  plane  of  incidence 
is  located  at  z  —  ct  =  -435  /mi.  The  line-outs  in  (c)  and  (d)  are  evaluated  at  y  =  0. 

CTR-like  feature  is  supposed  to  retain  the  form  of  the  vacuum  bunch  fields.  However,  as  is 
well  known,  dispersive  effects  lead  to  distortions  of  the  type  shown  in  the  figure. 

When  the  bunch  charge  is  high  enough  to  excite  large  nonlinear  polarization  currents  in 
the  material,  the  fields  induced  in  the  crystal  may  be  further  distorted.  Such  a  case  is  shown 
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FIG.  8:  Conservation  of  energy  during  simulated  bunch-crystal  interaction. 


in  Figs.  7(b)  and  (d).  For  the  chosen  crystal  orientation,  the  nonlinear  polarization  is 

P(")  —  Xl23-^yex  +  ^Xl23^xEyey  (34) 


Hence,  a  y-polarized  THz  field  induces  an  x-polarized  field  with  second  harmonic  com¬ 
ponents.  This  field  then  mixes  with  the  original  y-polarized  field  to  produce  a  distorted 
y-polarized  field.  This  distortion  can  be  clearly  seen  by  comparing  panels  (c)  and  (d)  of 
Fig.  7.  Qualitatively,  the  effect  is  to  produce  narrower  features  in  the  waveform.  In  an 
electro-optic  sensing  application,  this  could  adversely  affect  the  fidelity  of  the  diagnostic. 

Finally,  consider  the  flow  of  energy  in  the  crystal  during  the  simulation.  For  the  specific 
form  of  the  dielectric  response  considered  in  this  report,  Poynting’s  theorem  can  be  written 
as 


Idn's+?s//r<£2+c2B2>+g 


dUp 

dt 


+  £v  —  o 


(35) 


where  the  first  term  is  the  energy  radiated  through  a  surface  dfl,  the  second  is  the  rate  of 
change  of  field  energy  in  the  volume  fl,  and  the  third  is  the  rate  of  change  of  the  energy  of 
all  the  effective  particles  in  fh  Here,  the  effective  particle  energy  is 


Up  —  ^ mP 


la)  +9-0-(r-) 


(36) 


where  mp  is  the  mass,  rp  is  the  displacement,  qp  is  the  charge,  and  4>p(rp )  is  the  microscopic 
electrostatic  potential  of  the  displaced  particle.  The  rate  of  frictional  damping  for  an  effective 


particle  is 
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Cp  =  2Tpmp  (^j  (37) 

Fig.  8  displays  the  three  terms  in  Poynting’s  theorem,  integrated  over  time,  along  with 
the  sum.  The  “radiated”  curve  is  negative,  corresponding  to  the  fact  that  the  bunch  fields 
propagate  into  the  crystal.  The  steep  initial  slope  corresponds  to  the  generation  of  the  CTR- 
like  feature,  while  the  more  gradual  slope  corresponds  to  the  generation  of  the  Cherenkov-like 
feature.  The  “held”  and  “material”  curves  show  that  half  the  incident  energy  goes  to  the 
fields,  while  the  other  half  goes  to  the  effective  particles.  The  “sum”  curve  shows  that  energy 
is  conserved,  to  the  extent  that  it  vanishes.  One  notices  a  slight  discrepancy  at  the  moment 
the  bunch  fields  first  enter  the  crystal.  This  may  be  related  to  the  numerical  process  of 
generating  a  reflected  wave. 


VII.  CONCLUSIONS 

The  physics  of  nonlinear  crystals  excited  by  laser  pulses,  plasmas,  particle  beams,  or  any 
combination  of  these,  can  be  numerically  simulated  using  an  effective  particle  model  based 
on  particle-in-cell  techniques.  In  this  model,  bound  charges  respond  to  a  superposition 
of  microscopic  and  macroscopic  fields,  where  the  former  are  chosen  to  satisfy  the  known 
dispersion  characteristics  of  the  material,  and  the  latter  are  self-consistently  computed  using 
the  usual  FDTD  technique.  The  resulting  numerical  model  is  extremely  flexible  and  makes 
very  few  assumptions.  It  is,  however,  computationally  expensive.  Nevertheless,  it  is  now 
possible  to  carry  out  three  dimensional  simulations  of  electro-optic  sensing,  where  the  bunch 
holds  induced  in  an  electro-optic  crystal  are  computed  with  unprecedented  detail.  Other 
applications  of  this  model  might  include  modeling  the  nonlinear  optics  of  nano-composites, 
or  frequency  mixing  with  ultra-short  pulses. 
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Appendix 


In  this  appendix,  the  relationship  between  the  nonlinear  refractive  index,  ri2,  and  the  sec¬ 
ond  anharmonic  coefficient,  b ,  is  derived.  For  simplicity,  all  quantities  are  reduced  to  scalars. 
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Consider  a  perturbation  expansion  of  the  anharmonic  oscillator  equation.  Suppressing  the 
species  index,  and  expressing  all  quantities  in  frequency  space,  the  first  order  equation  is 

X>(w)r(1)  =  —E(u)  (38) 

m 

and  the  third  order  equation  is 

V(u j)r^  =  b  *  r^)  *  (39) 

Here,  the  asterisk  denotes  convolution  over  frequency,  and  the  electric  held  is  in  the  form 

(40) 

Inserting  r<'v>  into  the  third  order  equation  gives 


E{v)  =  -y  [£(w  -  W0)  +  6(w  +  w0)] 


d3)  =  & 


qE0\  3  F(u) 
2 m  )  V(l u) 


where 


_  —  3u>o)  3 5{lo  —  Wo)  3h(cu  +  Wo)  5{u  +  3cuo) 

JJ  )  ✓T'V  O  /  \  I  ✓T'V  O  /  \  ✓f'v  /  \  I  ✓T'*  O  /  \  ✓T'*  /  \  H- 


^3(w0)  D2(wo)P(-w0)  P2(-w0)P(w0)  7>3(-w0) 

Restoring  species  indices,  the  nonlinear  polarization  is 


P(3)  =  ^2N8q8r 
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This  quantity  must  now  be  connected  to  the  refractive  index. 

The  refractive  index  is  defined  by  n  =  e1/2,  with  e  the  relative  permittivity.  Using 
D  =  6q eE,  this  can  be  written  as 


n  = 


e0E  +  P 
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For  a  weakly  nonlinear  interaction,  P ^  -C  P and 

p{ 3) 

n  «  n0  + 
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where  no  is  the  linear  index,  given  by  Uq  =  1  +  P^/e^E.  The  nonlinear  index,  n2 ,  is 
conventionally  defined  by 

n  =  n0  +  n2/o  (46) 

where  Jo  is  the  cycle  averaged  intensity.  Choosing  the  Fourier  transform  convention  such 

that  Eq  is  the  peak  value  of  the  electric  held,  one  has  the  relationship  Jo  =  no£,(2/2r/o. 

Finally,  equating  Eqs.  (45)  and  (46),  and  evaluating  at  w0,  gives 

/  x  =  3  77o  y-  q2s _ WpS  , 
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